d <-#read_csv(paste0(home, "/data/for-analysis/dat.csv")) |> read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for-analysis/dat.csv") |>filter(age_ring =="Y") # use only length-at-age by filtering on age_ring#> Rows: 364546 Columns: 11#> ── Column specification ────────────────────────────────────────────────────────#> Delimiter: ","#> chr (5): age_ring, area, gear, ID, sex#> dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length#> #> ℹ Use `spec()` to retrieve the full column specification for this data.#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#> filter: removed 28,878 rows (8%), 335,668 rows remaining# Sample size by area and cohortns <- d |>group_by(cohort, area) |>summarise(n =n())#> group_by: 2 grouping variables (cohort, area)#> summarise: now 505 rows and 3 columns, one group variable remaining (cohort)# Minimum number of observations per area and cohortd$area_cohort <-as.factor(paste(d$area, d$cohort))d <- d |>group_by(area_cohort) |>filter(n() >100)#> group_by: one grouping variable (area_cohort)#> filter (grouped): removed 2,844 rows (1%), 332,824 rows remaining# Minimum number of observations per area, cohort, and aged$area_cohort_age <-as.factor(paste(d$area, d$cohort, d$age_bc))d <- d |>group_by(area_cohort_age) |>filter(n() >10)#> group_by: one grouping variable (area_cohort_age)#> filter (grouped): removed 3,974 rows (1%), 328,850 rows remaining# Minimum number of cohorts in a given areacnt <- d |>group_by(area) |>summarise(n =n_distinct(cohort)) |>filter(n >=10)#> group_by: one grouping variable (area)#> summarise: now 12 rows and 2 columns, ungrouped#> filter: no rows removedd <- d[d$area %in% cnt$area, ]# Plot cleaned dataggplot(d, aes(age_bc, length_mm)) +geom_jitter(size =0.1, height =0, alpha =0.1) +scale_x_continuous(breaks =seq(20)) +theme(axis.text.x =element_text(angle =0)) +theme(axis.text =element_text(size =12), axis.title =element_text(size =15)) +labs(x ="Age", y ="Length (mm)") +guides(color ="none") +facet_wrap(~area, scale ="free_x")
Code
# Longitude and latitude attributes for each areaarea <-c("BS", "BT", "FB", "FM", "HO", "JM", "MU", "RA", "SI_EK", "SI_HA", "TH", "VN")nareas <-length(area)lat <-c(60, 60.4, 60.3, 60.5, 63.7, 58, 59, 65.9, 57.3, 57.4, 56.1, 57.5)lon <-c(21.5, 18.1, 19.5, 18, 20.9, 16.8, 18.1, 22.3, 16.6, 16.7, 15.9, 16.9)area_attr <-data.frame(cbind(area = area, lat = lat, lon = lon)) |>mutate_at(c("lat","lon"), as.numeric)#> mutate_at: converted 'lat' from character to double (0 new NA)#> converted 'lon' from character to double (0 new NA)
Fit VBGE models
Code
# Get individual growth parameters (functions: VBGF/Gompertz,nls_out,fit_nls)IVBG <- d |>group_by(ID) |>summarize(k =nls_out(fit_nls(length_mm, age_bc, min_nage =4, model ="VBGF"))$k,k_se =nls_out(fit_nls(length_mm, age_bc, min_nage =4, model ="VBGF"))$k_se,linf =nls_out(fit_nls(length_mm, age_bc, min_nage =4, model ="VBGF"))$linf,linf_se =nls_out(fit_nls(length_mm, age_bc, min_nage =4, model ="VBGF"))$linf_se)#> group_by: one grouping variable (ID)#> summarize: now 91,422 rows and 5 columns, ungrouped
Inspect predictions
Code
IVBG <- IVBG |>drop_na(k) # The NAs are because the model didn't converge or because they were below the threshold age#> drop_na: removed 50,514 rows (55%), 40,908 rows remainingIVBG <- IVBG |>mutate(k_lwr = k -1.96*k_se,k_upr = k +1.96*k_se,linf_lwr = linf -1.96*linf_se,linf_upr = linf +1.96*linf_se,row_id =row_number())#> mutate: new variable 'k_lwr' (double) with 40,864 unique values and 0% NA#> new variable 'k_upr' (double) with 40,864 unique values and 0% NA#> new variable 'linf_lwr' (double) with 40,864 unique values and 0% NA#> new variable 'linf_upr' (double) with 40,864 unique values and 0% NA#> new variable 'row_id' (integer) with 40,908 unique values and 0% NA# Plot all K'sIVBG |>#filter(row_id < 5000) |>ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr)) +geom_point(alpha =0.5) +geom_errorbar(alpha =0.5) +NULL
# We can also consider removing the upper 95% quantile of standard errors (not quantile of K directly)IVBG |> tidylog::mutate(keep =ifelse(k >quantile(k_se, probs =0.95), "N", "Y")) |>#filter(row_id < 10000) |>ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr, color = keep)) +geom_point(alpha =0.5) +facet_wrap(~keep, ncol =1) +geom_errorbar(alpha =0.5) +NULL#> mutate: new variable 'keep' (character) with 2 unique values and 0% NA
Code
# Add back cohort and area variablesIVBG <- IVBG |>left_join(d |>select(ID, area_cohort)) |>separate(area_cohort, into =c("area", "cohort"), remove =TRUE, sep =" ") |>mutate(cohort =as.numeric(cohort))#> Adding missing grouping variables: `area_cohort_age`#> select: dropped 10 variables (length_mm, age_bc, age_catch, age_ring, area, …)#> Joining with `by = join_by(ID)`#> left_join: added 2 columns (area_cohort_age, area_cohort)#> > rows only in x 0#> > rows only in y (115,512)#> > matched rows 213,338 (includes duplicates)#> > =========#> > rows total 213,338#> mutate: converted 'cohort' from character to double (0 new NA)# Summarise and save for sample sizesample_size <- IVBG |>group_by(area) |>summarise(n_cohort =length(unique(cohort)),min_cohort =min(cohort),max_cohort =min(cohort),n_individuals =length(unique(ID)),n_data_points =n())#> group_by: one grouping variable (area)#> summarise: now 12 rows and 6 columns, ungroupedsample_size#> # A tibble: 12 × 6#> area n_cohort min_cohort max_cohort n_individuals n_data_points#> <chr> <int> <dbl> <dbl> <int> <int>#> 1 BS 18 1985 1985 2370 11832#> 2 BT 38 1972 1972 2193 10515#> 3 FB 48 1969 1969 8678 47933#> 4 FM 54 1963 1963 8660 48242#> 5 HO 34 1982 1982 2711 12454#> 6 JM 63 1953 1953 7738 40229#> 7 MU 19 1981 1981 1920 9930#> 8 RA 15 1985 1985 1205 5654#> 9 SI_EK 35 1968 1968 1396 6523#> 10 SI_HA 41 1967 1967 2991 15143#> 11 TH 5 2000 2000 85 340#> 12 VN 12 1987 1987 961 4543sample_size |>ungroup() |>summarise(sum_ind =sum(n_individuals), sum_n =sum(n_data_points))#> ungroup: no grouping variables#> summarise: now one row and 2 columns, ungrouped#> # A tibble: 1 × 2#> sum_ind sum_n#> <int> <int>#> 1 40908 213338write_csv(sample_size, paste0(home, "/output/sample_size.csv"))# Compare how the means and quantiles differ depending on this filteringIVBG_filter <- IVBG |> tidylog::drop_na(k_se) |> tidylog::filter(k_se <quantile(k_se, probs =0.95))#> drop_na: no rows removed#> filter: removed 10,670 rows (5%), 202,668 rows remaining# Summarize growth coefficients by cohort and areaVBG <- IVBG |>group_by(cohort, area) |>summarize(k_mean =mean(k, na.rm = T),k_median =quantile(k, prob =0.5, na.rm = T),linf_median =quantile(linf, prob =0.5, na.rm = T),k_lwr =quantile(k, prob =0.05, na.rm = T),k_upr =quantile(k, prob =0.95, na.rm = T)) |>ungroup()#> group_by: 2 grouping variables (cohort, area)#> summarize: now 382 rows and 7 columns, one group variable remaining (cohort)#> ungroup: no grouping variablesVBG_filter <- IVBG_filter |>group_by(cohort, area) |>summarize(k_mean =mean(k, na.rm = T),k_median =quantile(k, prob =0.5, na.rm = T),k_lwr =quantile(k, prob =0.05, na.rm = T),k_upr =quantile(k, prob =0.95, na.rm = T)) |>ungroup()#> group_by: 2 grouping variables (cohort, area)#> summarize: now 382 rows and 6 columns, one group variable remaining (cohort)#> ungroup: no grouping variablesggplot() +geom_ribbon(data = VBG, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr,fill ="All k's"), alpha =0.4) +geom_ribbon(data = VBG_filter, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr,fill ="Filtered"), alpha =0.4) +geom_line(data = VBG, aes(cohort, k_median, color ="All k's")) +geom_line(data = VBG_filter, aes(cohort, k_median, color ="Filtered")) +guides(fill =FALSE) +facet_wrap(~area)#> Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as#> of ggplot2 3.3.4.
Code
ggplot() +geom_line(data = VBG, aes(cohort, k_median, color ="All k's")) +geom_line(data = VBG_filter, aes(cohort, k_median, color ="Filtered")) +guides(fill =FALSE) +facet_wrap(~area)
Code
# No difference at all really. We should probably just discuss that with this model, achieving biologically reasonable parameter values and a good fit to data are sometimes two different things. In our case, we just want a representative value of the growth (as in time to reach average maximum size in the population) that accounts for the entire growth trajectory of an individual, for each area and cohort.
Add GAM-predicted temperature to growth models
Code
pred_temp <-read_csv(paste0(home, "/output/gam_predicted_temps.csv")) |>rename(cohort = year)#> Rows: 566 Columns: 5#> ── Column specification ────────────────────────────────────────────────────────#> Delimiter: ","#> chr (2): area, type#> dbl (3): year, temp, temp_gs#> #> ℹ Use `spec()` to retrieve the full column specification for this data.#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#> rename: renamed one variable (cohort)VBG <- VBG |>left_join(pred_temp, by =c("area", "cohort"))#> left_join: added 3 columns (temp, type, temp_gs)#> > rows only in x 0#> > rows only in y (184)#> > matched rows 382#> > =====#> > rows total 382# save data for map-plotcohort_sample_size <- IVBG |>group_by(area, cohort) |>summarise(n =n()) # individuals, not samples!#> group_by: 2 grouping variables (area, cohort)#> summarise: now 382 rows and 3 columns, one group variable remaining (area)VBG <-left_join(VBG, cohort_sample_size, by =c("cohort", "area"))#> left_join: added one column (n)#> > rows only in x 0#> > rows only in y ( 0)#> > matched rows 382#> > =====#> > rows total 382write_csv(VBG, paste0(home, "/output/vbg.csv"))
Plot VBGE fits
Code
# Sample 50 IDs per area and plot their data and VBGE fitsset.seed(4)ids <- IVBG |>distinct(ID, .keep_all =TRUE) |>group_by(area) |>slice_sample(n =30)#> distinct: removed 172,430 rows (81%), 40,908 rows remaining#> group_by: one grouping variable (area)#> slice_sample (grouped): removed 40,548 rows (99%), 360 rows remaining#ids |> ungroup() |> group_by(area) |> summarise(n = length(unique(ID))) |> arrange(n)IVBG2 <- IVBG |>filter(ID %in% ids$ID) |>distinct(ID, .keep_all =TRUE) |>select(ID, k, linf)#> filter: removed 211,555 rows (99%), 1,783 rows remaining#> distinct: removed 1,423 rows (80%), 360 rows remaining#> select: dropped 10 variables (k_se, linf_se, k_lwr, k_upr, linf_lwr, …)d2 <- d |>ungroup() |>filter(ID %in% ids$ID) |>left_join(IVBG2, by ="ID") |>mutate(length_mm_pred = linf*(1-exp(-k*age_bc)))#> ungroup: no grouping variables#> filter: removed 327,067 rows (99%), 1,783 rows remaining#> left_join: added 2 columns (k, linf)#> > rows only in x 0#> > rows only in y ( 0)#> > matched rows 1,783#> > =======#> > rows total 1,783#> mutate: new variable 'length_mm_pred' (double) with 1,783 unique values and 0% NAggplot(d2, aes(age_bc, length_mm, group = ID, color = ID)) +geom_jitter(width =0.3, height =0, alpha =0.5, size =0.4) +geom_line(aes(age_bc, length_mm_pred, group = ID, color = ID),inherit.aes =FALSE, alpha =0.8, linewidth =0.3) +guides(color ="none") +scale_color_viridis(discrete =TRUE, name ="Area", option ="cividis") +labs(x ="Age (yrs)", y ="Length (mm)") +facet_wrap(~factor(area, levels = (arrange(area_attr, desc(lat)))$area), ncol =3)
Code
ggsave(paste0(home, "/figures/vb_fits.pdf" ), width =17, height =17, units ="cm")rugs <- IVBG |>group_by(area) |>summarise(median_k =median(k),median_linf =median(linf))#> group_by: one grouping variable (area)#> summarise: now 12 rows and 3 columns, ungroupedk <- IVBG |>ggplot(aes(k, color =factor(area, order$area))) +geom_density(alpha =0.4, fill =NA, adjust =1.5) +scale_color_manual(values = colors, name ="Area") +coord_cartesian(expand =0) +labs(x =expression(italic(k))) +geom_rug(data = rugs, aes(median_k)) +theme(legend.position =c(0.9, 0.5))linf <- IVBG |>ggplot(aes(linf, color =factor(area, order$area))) +geom_density(alpha =0.4, fill =NA, adjust =1.5) +scale_color_manual(values = colors, name ="Area") +coord_cartesian(expand =0, xlim =c(0, 2000)) +labs(x =expression(paste(italic(L[infinity]), " [cm]"))) +geom_rug(data = rugs, aes(median_linf)) +guides(color ="none")k / linf
Code
ggsave(paste0(home, "/figures/vb_pars.pdf" ), width =17, height =17, units ="cm")
Fit Sharpe-Schoolfield model to K
By area
Code
model <-'sharpeschoolhigh_1981'# Get starting values on full dataset for Sharpe-Schoolfield modeldat <- VBG |>select(k_median, temp, area) |>rename(rate = k_median)#> select: dropped 8 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)lower <-get_lower_lims(dat$temp, dat$rate, model_name = model)upper <-get_upper_lims(dat$temp, dat$rate, model_name = model)start <-get_start_vals(dat$temp, dat$rate, model_name = model)# Sharpe-Schoolfield model fit to data for each areapreds <-NULLpred <-NULLfor(a inunique(dat$area)) {# Get data dd <- dat |>filter(area == a)# Fit model fit <-nls_multstart( rate ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dd,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new data new_data <-data.frame(temp =seq(min(dd$temp), max(dd$temp), length.out =100)) pred <-augment(fit, newdata = new_data) |>mutate(area = a)# Add to general data frame preds <-data.frame(rbind(preds, pred))}#> filter: removed 319 rows (84%), 63 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 328 rows (86%), 54 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 341 rows (89%), 41 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 347 rows (91%), 35 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 334 rows (87%), 48 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 344 rows (90%), 38 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 363 rows (95%), 19 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 348 rows (91%), 34 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 364 rows (95%), 18 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 367 rows (96%), 15 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 370 rows (97%), 12 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 377 rows (99%), 5 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA
All areas pooled
Code
# One sharpe schoolfield fitted to all datafit_all <-nls_multstart( rate ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dat,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new datanew_data_all <-data.frame(temp =seq(min(dat$temp), max(dat$temp), length.out =100))pred_all <-augment(fit_all, newdata = new_data_all) |>mutate(area ="all")#> mutate: new variable 'area' (character) with one unique value and 0% NA# Add t_optkb <-8.62e-05data.frame(par =names(coef(fit_all)), est =coef(fit_all)) |>pivot_wider(names_from = par, values_from = est) |>summarise(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1)))#> pivot_wider: reorganized (par, est) into (r_tref, e, eh, th) [was 4x2, now 1x4]#> summarise: now one row and one column, ungrouped#> # A tibble: 1 × 1#> t_opt#> <dbl>#> 1 9.72
Plot Sharpe Schoolfield fits
Code
# Plot growth coefficients by year and area against mean SSTp1 <- preds |>ggplot(aes(temp, .fitted, color =factor(area, order$area))) +geom_point(data = dat, aes(temp, rate, color =factor(area, order$area)), size =0.6) +geom_line(linewidth =1) +geom_line(data = pred_all, aes(temp, .fitted), linewidth =1, inherit.aes =FALSE, linetype =2) +scale_color_manual(values = colors, name ="Area") +guides(color =guide_legend(override.aes =list(size =1))) +scale_x_continuous(breaks =seq(-5, 20, 1)) +labs(x ="Temperature (°C)",y ="von Bertalanffy growth coefficient (*k*)") +theme(axis.title.y = ggtext::element_markdown())p1 +facet_wrap(~factor(area, levels = (arrange(area_attr, desc(lat)))$area))
Code
p1
Code
ggsave(paste0(home, "/figures/sharpe_school.pdf" ), width =15, height =15, units ="cm")
Extra analysis
Temperature in growing season instead of all year average
Fit Sharpe-Schoolfield model to K
By area
Code
model <-'sharpeschoolhigh_1981'# Get starting values on full dataset for Sharpe-Schoolfield modeldat2 <- VBG |>select(k_median, temp_gs, area) |>rename(rate = k_median,temp = temp_gs) # growing season! not annual average!#> select: dropped 8 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed 2 variables (rate, temp)lower <-get_lower_lims(dat2$temp, dat2$rate, model_name = model)upper <-get_upper_lims(dat2$temp, dat2$rate, model_name = model)start <-get_start_vals(dat2$temp, dat2$rate, model_name = model)# Sharpe-Schoolfield model fit to data for each areapreds2 <-NULLpred2 <-NULLfor(a inunique(dat2$area)) {# Get data dd <- dat2 |>filter(area == a)# Fit model fit2 <-nls_multstart( rate ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dd,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new data new_data2 <-data.frame(temp =seq(min(dd$temp), max(dd$temp), length.out =100)) pred2 <-augment(fit2, newdata = new_data2) |>mutate(area = a)# Add to general data frame preds2 <-data.frame(rbind(preds2, pred2))}#> filter: removed 319 rows (84%), 63 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 328 rows (86%), 54 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 341 rows (89%), 41 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 347 rows (91%), 35 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 334 rows (87%), 48 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 344 rows (90%), 38 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 363 rows (95%), 19 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 348 rows (91%), 34 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 364 rows (95%), 18 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 367 rows (96%), 15 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 370 rows (97%), 12 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 377 rows (99%), 5 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA
All areas pooled
Code
# One sharpe schoolfield fitted to all datafit_all2 <-nls_multstart( rate ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dat2,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new datanew_data_all2 <-data.frame(temp =seq(min(dat2$temp), max(dat2$temp), length.out =100))pred_all2 <-augment(fit_all2, newdata = new_data_all2) |>mutate(area ="all")#> mutate: new variable 'area' (character) with one unique value and 0% NA# Add t_optkb <-8.62e-05data.frame(par =names(coef(fit_all2)), est =coef(fit_all2)) |>pivot_wider(names_from = par, values_from = est) |>summarise(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1)))#> pivot_wider: reorganized (par, est) into (r_tref, e, eh, th) [was 4x2, now 1x4]#> summarise: now one row and one column, ungrouped#> # A tibble: 1 × 1#> t_opt#> <dbl>#> 1 22.1
Plot Sharpe Schoolfield fits
Code
# Plot growth coefficients by year and area against mean SSTpreds2 |>ggplot(aes(temp, .fitted, color =factor(area, order$area))) +geom_point(data = dat2, aes(temp, rate, color =factor(area, order$area)), size =0.6) +geom_line(linewidth =1) +geom_line(data = pred_all2, aes(temp, .fitted), linewidth =1, inherit.aes =FALSE, linetype =2) +scale_color_manual(values = colors, name ="Area") +guides(color =guide_legend(override.aes =list(size =1))) +scale_x_continuous(breaks =seq(-5, 20, 1)) +labs(x ="Temperature (°C)",y ="von Bertalanffy growth coefficient (*k*)") +theme(axis.title.y = ggtext::element_markdown())
Linf instead of k
By area
Code
model <-'sharpeschoolhigh_1981'# Get starting values on full dataset for Sharpe-Schoolfield modeldat3 <- VBG |>select(linf_median, temp, area) |>rename(rate = linf_median)#> select: dropped 8 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)lower <-get_lower_lims(dat3$temp, dat3$rate, model_name = model)upper <-get_upper_lims(dat3$temp, dat3$rate, model_name = model)start <-get_start_vals(dat3$temp, dat3$rate, model_name = model)# Sharpe-Schoolfield model fit to data for each areapreds3 <-NULLpred3 <-NULLfor(a inunique(dat3$area)) {# Get data dd <- dat3 |>filter(area == a)# Fit model fit3 <-nls_multstart( rate ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dd,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new data new_data3 <-data.frame(temp =seq(min(dd$temp), max(dd$temp), length.out =100)) pred3 <-augment(fit3, newdata = new_data3) |>mutate(area = a)# Add to general data frame preds3 <-data.frame(rbind(preds3, pred3))}#> filter: removed 319 rows (84%), 63 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 328 rows (86%), 54 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 341 rows (89%), 41 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 347 rows (91%), 35 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 334 rows (87%), 48 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 344 rows (90%), 38 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 363 rows (95%), 19 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 348 rows (91%), 34 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 364 rows (95%), 18 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 367 rows (96%), 15 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 370 rows (97%), 12 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 377 rows (99%), 5 rows remaining#> mutate: new variable 'area' (character) with one unique value and 0% NA
All areas pooled
Code
# One sharpe schoolfield fitted to all datafit_all3 <-nls_multstart( rate ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dat3,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new datanew_data_all3 <-data.frame(temp =seq(min(dat3$temp), max(dat3$temp), length.out =100))pred_all3 <-augment(fit_all3, newdata = new_data_all3) |>mutate(area ="all")#> mutate: new variable 'area' (character) with one unique value and 0% NA
Plot Sharpe Schoolfield fits
Code
# Plot growth coefficients by year and area against mean SSTpreds3 |>ggplot(aes(temp, .fitted, color =factor(area, order$area))) +geom_point(data = dat3, aes(temp, rate, color =factor(area, order$area)), size =0.6) +geom_line(linewidth =1) +geom_line(data = pred_all3, aes(temp, .fitted), linewidth =1, inherit.aes =FALSE, linetype =2) +scale_color_manual(values = colors, name ="Area") +guides(color =guide_legend(override.aes =list(size =1))) +scale_x_continuous(breaks =seq(-5, 20, 1)) +labs(x ="Temperature (°C)",y =expression(paste(italic(L[infinity]), " [cm]"))) +theme(axis.title.y = ggtext::element_markdown())
